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Abstract  Extinction  of  an  epidemic  or  a  species  is  a  rare  event  that  occurs  due  to  a  large, 
rare  stochastic  fluctuation.  Although  the  extinction  process  is  dynamically  unstable,  it  fol¬ 
lows  an  optimal  path  that  maximizes  the  probability  of  extinction.  We  show  that  the  opti¬ 
mal  path  is  also  directly  related  to  the  finite-time  Lyapunov  exponents  of  the  underlying 
dynamical  system  in  that  the  optimal  path  displays  maximum  sensitivity  to  initial  condi¬ 
tions.  We  consider  several  stochastic  epidemic  models,  and  examine  the  extinction  process 
in  a  dynamical  systems  framework.  Using  the  dynamics  of  the  finite-time  Lyapunov  ex¬ 
ponents  as  a  constructive  tool,  we  demonstrate  that  the  dynamical  systems  viewpoint  of 
extinction  evolves  naturally  toward  the  optimal  path. 

Keywords  Stochastic  dynamical  systems  and  Lyapunov  exponents  •  Optimal  path  to 
extinction 


1.  Introduction 

Control  and  eradication  of  infectious  diseases  are  among  the  most  important  goals  for 
improving  public  health.  Although  the  global  eradication  of  a  disease  (e.g.,  smallpox)  has 
been  achieved  (B reman  and  Arita,  1980),  it  is  difficult  to  accomplish  and  remains  a  public 
health  target  for  polio,  malaria,  and  many  other  diseases,  including  childhood  diseases 
(Aylward  et  al.,  2000).  The  global  eradication  of  a  disease  is  not  the  only  type  of  disease 
extinction  process.  Lor  example,  a  disease  may  “fade  out”  or  become  locally  extinct  in  a 
region.  Since  the  extinction  is  local,  the  disease  may  be  reintroduced  from  other  regions 
(Grassly  et  al.,  2005).  Additionally,  extinction  may  also  occur  to  individual  strains  of  a 
multistrain  disease  (Minayev  and  Lerguson,  2009),  such  as  influenza  or  dengue  fever.  It 
is  worth  noting  that  the  extinction  process  is  also  of  interest  in  the  fields  of  evolution  and 
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ecology.  As  an  example,  bio-diversity  arises  from  the  interplay  between  the  introduction 
and  extinction  of  species  (Azaele  et  al.,  2006;  Banavar  and  Maritan,  2009). 

In  order  to  predict  the  dynamics  of  disease  outbreak  and  spread  as  well  as  to  im¬ 
plement  control  strategies  which  promote  disease  extinction,  one  must  investigate  how 
model  parameters  affect  the  probability  of  extinction.  For  example,  Dykman  et  al.  (2008) 
showed  that  disease  control  and  extinction  depend  on  both  epidemiological  and  sociolog¬ 
ical  parameters  determining  epidemic  growth  rate.  Additionally,  since  extinction  occurs 
in  finite  populations,  another  factor  in  determining  extinction  risk  is  the  local  commu¬ 
nity  size  (Bartlett,  1957,  1960;  Keeling  and  Grenfell,  1997;  Conlan  and  Grenfell,  2007). 
Further  complications  arise  from  the  fact  that  in  general,  stochastic  effects  cause  random 
transitions  within  the  discrete,  finite  populations.  These  stochastic  effects  may  be  intrin¬ 
sic  to  the  system  or  may  arise  from  the  external  environment.  Small  population  size,  low 
contact  frequency  for  frequency-dependent  transmission,  competition  for  resources,  and 
evolutionary  pressure  (de  Castro  and  Bolker,  2005),  as  well  as  heterogeneity  in  popula¬ 
tions  and  transmission  (Lloyd  et  al.,  2007)  may  all  be  determining  factors  for  extinction  to 
occur.  Other  factors  which  can  affect  the  risk  of  extinction  include  the  nature  and  strength 
of  the  noise  (Melbourne  and  Hastings,  2008),  disease  outbreak  amplitude  (Alonso  et  al., 
2006),  and  seasonal  phase  occurrence  (Stone  et  al.,  2007). 

In  large  populations,  the  intensity  of  the  intrinsic  noise  is  often  quite  small.  However, 
it  is  possible  that  a  rare,  large  fluctuation  which  occurs  with  finite  probability  can  cause 
the  system  to  reach  the  extinct  state.  Since  the  extinct  state  is  absorbing  due  to  effective 
stochastic  forces,  eventual  extinction  is  guaranteed  when  there  is  no  source  of  reintro¬ 
duction  (Bartlett,  1949;  Allen  and  Burgin,  2000;  Gardiner,  2004).  However,  because  fade 
outs  are  usually  rare  events  in  large  populations,  typical  time  scales  for  extinction  may  be 
extremely  long. 

Birth-death  systems  (Gardiner,  2004;  van  Kampen,  2007),  which  possess  intrinsic 
noise,  form  an  important  class  of  stochastic  processes.  These  systems  have  been  used 
in  the  field  of  mathematical  epidemiology  (Bartlett,  1961;  Andersson  and  Britton,  2000; 
Choisy  et  al.,  2007).  In  these  systems,  the  intrinsic  noise  arises  from  the  discreteness  of 
the  individuals  (or  species)  and  the  randomness  of  their  interactions.  To  predict  probabili¬ 
ties  of  events  occurring  at  certain  times,  a  description  of  the  stochastic  system  is  provided 
using  the  master  equation  formalism.  Stochastic  master  equations  are  commonly  used  in 
statistical  physics  when  dealing  with  chemical  reaction  processes  (Kubo,  1963). 

For  systems  with  a  large  number  of  individuals,  a  van  Kampen  system- size  expansion 
may  be  used  to  approximate  the  master  equation  by  a  Fokker-Planck  equation  (Gardiner, 
2004;  van  Kampen,  2007).  However,  the  technique  fails  in  determining  the  probability  of 
large  fluctuations  (Gaveau  et  al.,  1996;  Elgart  and  Kamenev,  2004).  Instead,  in  this  article, 
we  will  employ  an  eikonal  approximation  to  recast  the  problem  in  terms  of  an  effective 
classical  Hamiltonian  system  (Kamenev  and  Meerson,  2008).  With  high  probability,  the 
intrinsic  noise  will  induce  extinction  of  the  disease  or  species  along  a  heteroclinic  trajec¬ 
tory  in  the  phase  space  of  the  classical  Hamiltonian  flow.  This  trajectory  is  known  as  the 
optimal  path  (to  extinction). 

It  is  highly  desirable  to  locate  the  optimal  path  since  the  extinction  rate  depends  on 
the  probability  to  traverse  this  path.  Additionally,  the  effect  of  a  control  strategy  on  the 
extinction  rate  can  be  determined  by  its  effect  on  the  optimal  path  (Dykman  et  al.,  2008). 
Through  the  use  of  the  eikonal  approximation,  we  consider  a  mechanistic  dynamical  sys¬ 
tems  problem  rather  than  the  original  stochastic  problem.  Unlike  other  methods,  this  ap¬ 
proach  enables  one  both  to  estimate  extinction  lifetimes  and  to  draw  conclusions  about 
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the  path  taken  to  extinction.  This  more  detailed  understanding  of  how  extinction  occurs 
may  lead  to  new  stochastic  control  strategies  (Dykman  et  al.,  2008). 

In  general,  the  optimal  path  to  extinction  is  an  unstable  dynamical  object.  Therefore, 
many  researchers  have  investigated  the  scaling  of  extinction  rates  with  respect  to  a  pa¬ 
rameter  near  a  bifurcation  point,  where  the  dynamics  are  slow  (Doering  et  al.,  2005; 
Kamenev  and  Meerson,  2008;  Kamenev  et  al.,  2008;  Dykman  et  al.,  2008).  The  analytical 
calculation  of  extinction  rates  far  from  a  bifurcation  point  is  much  more  difficult,  and  thus 
researchers  often  resort  to  averaging  over  many  stochastic  runs  (e.g.,  Shaw  and  Schwartz, 
2010).  The  numerical  computation  of  the  optimal  path  trajectory  has  been  achieved  in  the 
past  using  a  shooting  method  (Kamenev  and  Meerson,  2008).  However,  since  the  proce¬ 
dure  is  very  sensitive  to  boundary  conditions,  it  is  difficult  to  implement  when  analyzing 
paths  far  away  from  bifurcation  points  or  when  analyzing  high-dimensional  models. 

In  this  article,  we  develop  a  novel  method  for  computing  the  optimal  extinction  tra¬ 
jectory  that  relies  on  the  calculation  of  the  dynamical  system’s  finite-time  Lyapunov  ex¬ 
ponents  (FTLE).  The  classical  Lyapunov  exponent  provides  a  quantitative  measure  of 
how  infinitesimally  close  particles  in  a  dynamical  system  behave  asymptotically  as  time 
t  ->  =boo  (Guckenheimer  and  Holmes,  1986).  Similarly,  the  FTLE  provides  a  quantita¬ 
tive  measure  of  how  much  nearby  particles  separate  after  a  specific  amount  of  time  has 
elapsed. 

The  FTLE  fields  were  used  by  Pierrehumbert  (1991)  and  Pierrehumbert  and  Yang 
(1993)  to  characterize  structures,  including  mixing  regions  and  transport  barriers,  in  the 
atmosphere.  Later,  these  structures  were  quantified  more  rigorously  by  introducing  the 
idea  of  Lagrangian  Coherent  Structures  (LCS)  (Haller,  2000,  2001,  2002;  Shadden  et  al., 
2005;  Lekien  et  al.,  2007;  Branicki  and  Wiggins,  2010).  The  definition  of  a  LCS  as  a 
ridge  of  the  FTLE  field  was  introduced  by  Haller  (2002)  and  formalized  by  Shadden  et  al. 
(2005).  The  FTLE  method  provides  a  measure  of  how  sensitively  the  system’s  future  be¬ 
havior  depends  on  its  current  state,  and  the  LCS,  or  ridge,  is  a  structure  which  has  locally 
maximal  FTLE  value.  In  this  article,  we  will  show  that  the  system  displays  maximum 
sensitivity  near  the  optimal  extinction  trajectory,  which  enables  us  to  dynamically  evolve 
toward  the  optimal  escape  trajectory  using  FTLE  calculations. 

In  this  article,  we  consider  three  standard  epidemic  models  that  contain  intrinsic  or 
external  noise  sources  and  illustrate  the  power  of  our  method  to  locate  the  optimal  path 
to  extinction.  Even  though  our  examples  are  taken  from  infectious  disease  models,  the 
approach  is  applicable  to  any  extinction  process  or  escape  process.  Section  2  discusses 
stochastic  modeling  in  the  limit  of  large  population  size,  while  Section  3  describes  the 
theory  that  underlies  the  Lyapunov  exponent  computations.  In  Section  4,  our  method  is 
used  to  find  the  optimal  path  to  extinction  for  three  examples.  In  Section  5,  we  demon¬ 
strate  that  the  optimal  path  to  extinction  also  possesses  a  local  maximum  of  the  FTLE 
field,  and  conclusions  are  contained  in  Section  6. 


2.  Stochastic  modeling  in  the  large  population  limit 

To  introduce  the  idea  of  constructing  an  optimal  path  in  stochastic  dynamical  systems,  we 
consider  the  problem  of  extinction  taken  from  mathematical  epidemiology.  The  stochastic 
formulation  of  the  problem  accounts  for  the  random  state  transitions  which  govern  the 
dynamics  of  the  system.  To  quantitatively  account  for  the  randomness  in  the  system,  we 
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will  formulate  a  master  equation  which  describes  the  time  evolution  of  the  probability  of 
finding  the  system  in  a  particular  state  at  a  certain  time  (Gardiner,  2004;  van  Kampen, 
2007). 

The  state  variables  X  g  W1  of  the  system  describe  the  components  of  a  population, 
while  the  random  state  transitions  which  govern  the  dynamics  are  described  by  the  tran¬ 
sition  rates  W(X,r),  where  r  e  W1  is  an  increment  in  the  change  of  X.  Consideration 
of  the  net  change  in  increments  of  the  state  of  the  system  results  in  the  following  master 
equation  for  the  probability  density  p(X,  t)  of  finding  the  system  in  state  X  at  time  t : 

9P(^’°  =  J2iW^X  ~  r:  r)P(X  ~r^~  W(V  r)p(X,  t)].  (1) 

r 

If  the  total  population  size  of  the  system  is  N,  the  components  of  the  state  variable 
can  be  rescaled  to  be  fractions  of  the  population  by  letting  x  =  X/N.  Then  the  general 
solution  (Kubo  et  al.,  1973)  for  the  transition  probability  from  x0  to  x  in  the  time  interval 
from  to  to  t  can  be  written  as  the  following  path  integral: 

P{x,  t\xo,  to)  =  J  dV(x,  p)  expj  —  N  <is[//(x(s),  p(s),  s)  —  /?(s)x(s)]  j,  (2) 

where  dV( x,  p)  denotes  integration  over  all  paths. 

The  integral  in  the  exponent  of  Eq.  (2)  is  the  action,  and  the  Hamiltonian  H(x,  p\  t)  is 
given  in  general  as 

H(x ,  p\  t )  =  w(x>  f*)(exp (p  •  r)  —  l),  (3) 


where  w(x\r)  is  defined  as  the  transition  rate  W  per  individual.  The  Hamiltonian  H 
depends  both  on  the  state  of  the  system  x  as  well  as  the  momentum  p,  which  provides  an 
effective  force  due  to  stochastic  fluctuations  on  the  system.  It  should  be  noted  that  instead 
of  using  the  Hamiltonian  representation,  one  could  use  the  Lagrangian  representation, 
which  results  in  the  following  alternative  solution: 


P(x,  t\xo,  to)  =  J  dV(x)exip\^N  ^  ds  L(x(s),  x(s),  s)  j, 


(4) 


with  L(x,  x;  t)  =  —H(x,  p\  t)  +  x  •  p. 

The  action  reveals  much  information  about  the  probability  evolution  of  the  system, 
from  scaling  near  bifurcation  points  in  non-Gaussian  processes  to  rates  of  extinction  as  a 
function  of  epidemiological  parameters  (Dykman,  1990;  Dykman  et  al.,  2008).  In  order 
to  maximize  the  probability  of  extinction,  one  must  minimize  the  action.  The  minimizing 
formulation  entails  finding  the  solution  to  the  Hamilton-Jacobi  equation,  which  means 
that  one  must  solve  the  2n -dimensional  system  of  Hamilton’s  equations  for  x  and  p. 
The  appropriate  boundary  conditions  of  the  system  are  such  that  a  solution  starts  at  a 
nonzero  state,  such  as  an  endemic  state,  and  asymptotically  approaches  one  or  more  zero 
components  of  the  state  vector.  Therefore,  a  trajectory  that  is  a  solution  to  the  two-point 
boundary  value  problem  determines  a  path,  which  in  turn  yields  the  probability  of  going 
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from  the  initial  state  to  the  final  state.  The  optimal  path  to  extinction  is  the  path  that 
minimizes  the  action  in  either  the  Hamiltonian  or  Lagrangian  representation. 

To  illustrate  an  application  of  the  theory  for  a  finite  population,  we  consider  in  de¬ 
tail  an  explicit  example,  namely  the  well-known  problem  of  extinction  in  a  Susceptible- 
Infectious-Susceptible  (SIS)  epidemic  model.  In  this  example,  the  population  consists  of 
S  susceptible  individuals  and  I  infectious  individuals.  The  population  is  driven  via  the 
birth  rate  /z,  which  is  also  equal  to  the  death  rate.  The  transition  rates  for  X  =  ( S ,  I)T  are 
given  as  follows: 

W(X\  (1,  0))  =  Nfi,  W(X\  (-1,0))  =  iiX i, 

W(X;  (0,  -1))  =  /iX2,  W(X;  (1,  -1))  =  yX2,  (5) 

W(X;  (-1,  1))=j8X1X2/N, 

where  /3  is  the  mass  action  contact  rate,  y  is  the  recovery  rate,  and  N  is  now  a  parameter 
for  the  average  size  of  the  population.  For  large  S,  I  oc  N,  fluctuations  of  S  and  I  are  small 
on  average.  If  these  fluctuations  are  disregarded,  one  arrives  at  the  following  deterministic 
mean-field  equations  for  S  and  I : 


Xl  =  N/i  —  iiX i  +  yX2  -  /3XiX2/N,  (6a) 

X2  =  -(/z  +  y)X2  +  PXxX2/N.  (6b) 

Equations  (6a)-(6b)  are  the  standard  equations  of  the  SIS  model  in  the  absence  of  fluc¬ 
tuations.  For  the  parameter  Rq  =  P/(fi  +  y)  >  1,  they  have  a  stable,  attracting  solution 
XA  =  Nxa  with  X\A  =  Rq1,  and  x2A  =  1  —  R q1,  which  corresponds  to  endemic  dis¬ 
ease.  In  addition,  Eqs.  (6a)-(6b)  have  an  unstable  stationary  state  (saddle  point)  given  by 
Xs  =  Nxs  with  x\s  =  1  and  x2 s  =  0,  which  corresponds  to  the  extinct,  or  disease-free, 
state. 

For  N  1,  the  steady  state  distribution  p(X)  has  a  peak  at  the  stable  state  XA  with 
width  a  Nl/1.  This  peak  is  formed  over  a  typical  relaxation  time  given  in  Dykman  et  al. 
(2008)  and  Schwartz  et  al.  (2009).  However,  in  the  process  of  extinction,  we  are  interested 
in  the  probability  of  having  a  small  number  of  infectious  individuals,  which  is  determined 
by  the  tail  of  the  distribution.  The  distribution  tail  can  be  obtained  by  seeking  the  solution 
of  Eqs.  (6a)-(6b)  in  the  eikonal  form  (Elgart  and  Kamenev,  2004;  Doering  et  al.,  2005; 
Kubo  et  al.,  1973;  Wentzell,  1976;  Gang,  1987;  Dykman  et  al.,  1994;  Tretiakov  et  al., 
2003)  given  by 

p{X)  =  exp[— Af<S(x)l,  x  =  X/N, 

(7) 

p{X  +  r)  «  p{X)  exp(— p  •  r),  p  =  dS(x)/dx, 
where  <S(x)  is  the  action. 

To  leading  order  in  N~l ,  the  equation  for  <S(x)  has  a  form  of  the  Hamilton- Jacobi 
equation  S  =  —H(x,  dxS;  t ),  where  S  is  the  effective  action,  and  the  effective  Hamil¬ 
tonian  is  given  by  Eq.  (3),  with  tz>(x;  r)  =  N~lW{X\  r )  being  the  transition  rates  per 
individual.  The  action  <S(x)  can  be  found  from  classical  trajectories  of  the  auxiliary  sys¬ 
tem  with  Hamiltonian  H  that  satisfy  the  following  equations: 


x  =  dpH(x,  p ), 


p  =  -dxH(x,  p). 


(8) 
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Since  the  maximum  of  the  probability  of  extinction  is  found  by  minimizing  the  action, 
we  compute  the  trajectory  satisfying  the  Hamiltonian  system  that  has  as  its  asymptotic 
limits  in  time  the  endemic  state  as  t  —  oo  and  the  extinct  state  as  t  — ►  +oo.  The  ac¬ 
tion  then  has  the  form  from  Eq.  (2)  (Wentzell,  1976;  Gang,  1987;  Dykman  et  al.,  1994; 
Tretiakov  et  al.,  2003): 


(9) 


In  Eq.  (9),  the  integral  is  calculated  for  a  Hamiltonian  trajectory  (jc(f),  pit))7  that  starts 
as  t  —>  —oo  at  x  xA,  p  —>  0,  and  arrives  as  t  oo  at  the  state  xS-  This  trajectory 
describes  the  most  probable  sequence  of  elementary  events  x  —>  x  +  r  that  brings  the 
system  to  xs. 

Several  authors  have  considered  how  extinction  rates  scale  with  respect  to  a  pa¬ 
rameter  near  bifurcation  points  (Doering  et  al.,  2005;  Kamenev  and  Meerson,  2008; 
Kamenev  et  al.,  2008;  Dykman  et  al.,  2008)  when  the  distance  to  the  bifurcation  point 
is  small  and  the  dynamics  is  very  slow.  For  an  epidemic  model,  this  means  that  the  repro¬ 
ductive  rate  of  infection  is  greater  than  but  very  close  to  one.  However,  most  real  diseases 
have  reproductive  rates  of  infection  greater  than  1.5  (Anderson  and  May,  1991),  which 
translates  into  faster  growth  rates  from  the  extinct  state.  In  general,  in  order  to  get  ana¬ 
lytic  scaling  results,  one  must  compute  the  optimal  path  using  either  the  Hamiltonian  or 
Lagrangian  equations  of  motion.  However,  far  from  bifurcation  points,  one  is  seldom  able 
to  perform  the  required  analysis  or  computation. 

Additionally,  the  computation  of  the  optimal  path  involves  the  use  of  a  numerical 
shooting  method  (Kamenev  and  Meerson,  2008).  The  Hamiltonian  or  Lagrangian  repre¬ 
sentation  of  an  n -dimensional  dynamical  system  lies  in  2n -dimensional  space.  Therefore, 
for  even  relatively  low-dimensional  dynamical  systems,  the  use  of  a  shooting  method 
to  find  the  optimal  path  to  extinction  can  be  quite  problematic.  In  the  next  section,  we 
demonstrate  how  to  evolve  naturally  to  the  optimal  path  using  a  dynamical  systems  ap¬ 
proach. 


3.  Finite- time  Lyapunov  exponents 

We  consider  a  velocity  field  v  :  M2n  x  I  M2n  given  by  the  Hamiltonian  field  in  Eq.  (8) 
which  is  defined  over  the  time  interval  /  =  and  the  following  system  of  equa¬ 

tions: 


Ht;ti,y0)  =  v(y(t;ti,y0),t), 
y(ti',ti,y0)  =  y0 , 


(10a) 

(10b) 


where  y  =  (x,  p)T  e  R2”,  yQ  e  R2”,  and  t  €  I. 

Such  a  continuous  time  dynamical  system  has  quantities,  known  as  Lyapunov  expo¬ 
nents,  which  are  associated  with  the  trajectory  of  the  system  in  an  infinite  time  limit.  The 
Lyapunov  exponents  measure  the  growth  rates  of  the  linearized  dynamics  about  the  tra¬ 
jectory.  To  find  the  finite- time  Lyapunov  exponents  (FTLE),  one  computes  the  Lyapunov 
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exponents  on  a  restricted  finite  time  interval.  For  each  initial  condition,  the  exponents 
provide  a  measure  of  its  sensitivity  to  small  perturbations.  Therefore,  the  FTLE  is  a  mea¬ 
sure  of  the  local  sensitivity  to  initial  data.  For  the  purpose  of  completeness,  we  briefly 
recapitulate  the  derivation  of  the  FTLE.  Details  regarding  the  derivation  along  with  the 
appropriate  smoothness  assumptions  can  be  found  in  Haller  (2000,  2001,  2002),  Shadden 
et  al.  (2005),  Lekien  et  al.  (2007),  and  Branicki  and  Wiggins  (2010). 

The  integration  of  Eqs.  (10a)-(10b)  from  the  initial  time  tt  to  the  final  time  tt  +  T 

t- -\-T 

yields  the  flow  map  4>tl.  which  is  defined  as  follows: 

<frli+T  ■  J’o  ^  <P‘t‘+T(yo)  =  y(ti  +  T ;  ti,  y0).  (11) 


Then  the  FTLE  can  be  defined  as 

a(y,ti,T)=-^\n^/x^(A),  (12) 

where  A.max(Z\)  is  the  maximum  eigenvalue  of  the  right  Cauchy-Green  deformation  tensor 
A,  which  is  given  as  follows: 


A(y,tifT)  = 


V  dy(t) 


^iyov))\ 

dy(t)  )’ 


(13) 


with  *  denoting  the  adjoint. 

For  a  given  y  e  at  an  initial  time  tt,  Eq.  (12)  gives  the  maximum  finite-time  Lya¬ 
punov  exponent  for  some  finite  integration  time  T  (forward  or  backward),  and  provides  a 
measure  of  the  sensitivity  of  a  trajectory  to  small  perturbations. 

The  FTLE  field  given  by  a(y,  T )  can  be  shown  to  exhibit  “ridges”  of  local  max¬ 
ima  in  phase  space.  The  ridges  of  the  field  indicate  the  location  of  attracting  (backward 
time  FTLE  field)  and  repelling  (forward  time  FTLE  field)  structures.  In  two-dimensional 
(2D)  space,  the  ridge  is  a  curve  which  locally  maximizes  the  FTLE  field  so  that  trans¬ 
verse  to  the  ridge  one  finds  the  FTLE  to  be  a  local  maximum.  What  is  remarkable  is 
that  the  FTLE  ridges  correspond  to  the  optimal  path  trajectories,  which  is  shown  heuris- 
tically  in  Section  5.  The  basic  idea  is  that  since  the  optimal  path  is  inherently  unsta¬ 
ble  and  observed  only  through  many  realizations  of  stochastic  experiments,  the  FTLE 
shows  that  locally,  the  path  is  also  the  most  sensitive  to  initial  data.  Figure  1  shows  a 
schematic  that  demonstrates  why  the  optimal  path  has  a  local  maximum  to  sensitivity. 
If  one  chooses  an  initial  point  on  either  side  of  the  path  near  the  endemic  state,  the  two 
trajectories  will  separate  exponentially  in  time.  This  is  due  to  the  fact  that  both  extinct  and 
endemic  states  are  unstable,  and  the  connecting  trajectory  defining  the  path  is  unstable  as 
well.  Any  initial  points  starting  near  the  optimal  path  will  leave  the  neighborhood  in  short 
time. 


4.  Finding  the  optimal  path  to  extinction  using  FTLE 

We  now  apply  our  theory  of  dynamical  sensitivity  to  the  problem  of  locating  optimal 
paths  to  extinction  for  several  examples.  We  consider  the  case  of  internal  fluctuations, 
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Fig.  1  Schematic  showing  the  path  from  the  endemic  state  (blue)  to  the  extinct  state  (red).  The  optimal 
path  leaves  the  endemic  point  along  an  unstable  manifold  and  connects  to  the  extinct  state  along  a  stable 
manifold.  The  magenta  and  green  dashed  lines  represent  trajectories  on  either  side  of  the  optimal  path.  The 
initial  starting  distance  between  trajectories  near  the  endemic  state  expands  exponentially  in  forward  time 
(shown  by  the  cyan  lines).  Locally,  this  demonstrates  that  the  finite-time  Lyapunov  measure  of  sensitivity 
with  respect  to  initial  data  is  maximal  along  the  optimal  path.  This  is  evident  in  the  ridges  observed  in  the 
evolution  of  the  exponents.  (Color  figure  online.) 


where  the  noise  is  not  known  a  priori,  as  well  as  the  case  of  external  noise,  where  the 
noise  is  specified.  In  each  case,  the  interaction  of  the  noise  and  state  of  the  systems  begin 
with  a  description  of  the  Hamiltonian  (or  Lagrangian)  to  describe  the  unstable  flow.  Then 
the  corresponding  equations  of  motion  are  used  to  compute  the  ridges  corresponding  to 
maximum  FTLE,  which  in  turn  correspond  to  the  optimal  extinction  paths  (Schwartz  et 
al.,  2010). 

4.1.  Example  1 — Extinction  in  a  branching-annihilation  process 

For  an  example  of  a  system  with  intrinsic  noise  fluctuations  which  has  an  analytical  solu¬ 
tion,  we  consider  extinction  in  the  stochastic  branching-annihilation  process 

a4>2A  and  2 A  L  0,  (14) 

where  X  and  /i  >  0  are  constant  reaction  rates  (Elgart  and  Kamenev,  2004;  Assaf  et  al., 
2008).  Equation  (14)  is  a  single  species  birth-death  process  and  can  be  thought  of  as  a 
simplified  form  of  the  Verhulst  logistic  model  for  population  growth  (Nasell,  2001). 

The  stochastic  process  given  by  Eq.  (14)  contains  intrinsic  noise  which  arises  from 
the  randomness  of  the  reactions  and  the  fact  that  the  population  consists  of  discrete  in¬ 
dividuals.  This  intrinsic  noise,  which  can  generate  a  rare  sequence  of  events  that  cause 
the  system  to  evolve  to  the  empty  state,  can  be  accounted  for  using  a  master  equation 
approach.  The  probability  Pn(t)  to  observe,  at  time  t,  n  individuals  is  governed  by  the 
following  master  equation: 

Pn  =  “[(«  +  2)(«  +  P)Pn+2  —  n(n  -  \)Pn\  +  X[(n  -  l)P„_i  -  nP„\ 


(15) 
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Using  the  formalism  of  Assaf  et  al.  (2008),  Eq.  (15)  is  recast  as  the  following  exact 
evolution  equation  for  G(p,t): 


dG 

dt 


[i 


d2G 


(16) 


where  G  is  a  probability  generating  function  given  by 


oo 

G(p,t)  =  YJPnPn(.t), 

n= 0 


(17) 


and  where  p  is  an  auxiliary  variable. 

We  substitute  the  eikonal  ansatz  G(p,t)  =  exp[— S(p,  t)],  where  S  is  the  action,  into 
Eq.  (16)  and  neglect  the  higher-order  d2S/dp2  term.  This  results  in  a  Hamilton- Jacobi 
equation  for  S(p,  t).  By  introducing  a  conjugate  coordinate  q  =  —dS/dp  and  by  shifting 
the  momentum  p  =  p  —  1,  then  one  arrives  at  the  following  Hamiltonian: 


H(q,p)=yHl+p)-^(2  +  p)qjqp.  (18) 

Hamilton’s  equations  are  therefore  given  as 

q  =  ^-=q[X(l+2p)-n(l+p)q],  (19a) 

p  = -^- =  p[/x(2  +  p)q  -  +  p)].  (19b) 

The  Hamiltonian  given  by  Eq.  (18)  has  three  zero-energy  curves.  The  first  is  the  mean- 
field  zero-energy  line  p  =  0,  which  contains  two  hyperbolic  points  given  as  h\  =  (q,  p)  = 
(A //x,  0)  and  ho  =  (q,  p)  =  (0,  0).  The  second  is  the  extinction  line  q  =  0,  which  contains 
another  hyperbolic  point  given  as  h2  =  (q,  p)  =  (0,-1).  The  third  zero-energy  curve  is 
non- trivial  and  is  given  as  follows: 


2A(1  +  p) 
/x(2  +  /?)  ' 


(20) 


The  segment  of  the  curve  given  by  Eq.  (20)  which  lies  between  —  i  <p<  0  corresponds 
to  a  heteroclinic  trajectory.  As  t  approaches  —  oo,  the  trajectory  exits  the  hyperbolic  point 
h  i  along  its  unstable  manifold  and  enters  the  hyperbolic  point  along  its  stable  mani¬ 
fold  as  t  approaches  oo.  This  heteroclinic  trajectory  is  the  optimal  path  to  extinction  and 
describes  the  most  probable  (rare)  sequence  of  events  which  evolves  the  system  from  a 
quasi- stationary  state  to  extinction  (Assaf  et  al.,  2008). 

To  show  that  the  FTLE  evolves  to  the  optimal  path,  we  calculate  the  FTLE  field  using 
Eqs.  (19a)-(19b).  Figure  2(a)  shows  the  forward  FTLE  plot  computed  using  Eqs.  (19a)- 
(19b)  for  T  =  6,  with  A  =  2.0  and  fi  =  0.5.  In  Fig.  2(a),  one  can  see  that  the  optimal  path 
to  extinction  is  given  by  the  ridge  associated  with  the  maximum  FTLE. 
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Fig.  2  FTLE  flow  fields  computed  using  Eqs.  (19a)-(19b)  with  X  =  2.0  and  /x  =  0.5.  The  integration  time 
is  T  =  6  with  an  integration  step  size  of  t  =  0.1  and  a  grid  resolution  of  0.01  in  both  q  and  p.  (a)  Forward 
FTLE  field  with  the  optimal  path  to  extinction  given  by  the  LCS.  (b)  Average  of  the  forward  and  backward 
FTLE  fields  with  the  three  zero-energy  curves  given  by  the  LCS  and  overlaid  with  the  analytical  solution 
of  these  curves  given  by  p  =  0,  q  =  0,  and  Eq.  (20).  Note  that  the  averaging  affects  only  the  value  of  the 
FTLE  and  not  the  structure  of  the  FTLE  field.  (Color  figure  online.) 


Not  all  of  the  attracting  structures  are  shown  in  Fig.  2(a)  because  the  maximum  for¬ 
ward  time  FTLE  identifies  repelling  structures  where  nearby  initial  conditions  diverge. 
The  other  attracting  structures  may  be  found  by  computing  the  backward  FTLE  field. 
By  overlaying  the  forward  and  backward  FTLE  fields,  one  can  see  in  their  entirety  all 
three  zero-energy  curves  including  the  optimal  path  to  extinction  in  Fig.  2(b).  Also  shown 
in  Fig.  2(b)  are  the  analytical  solutions  to  the  three  zero-energy  curves  given  by  p  =  0, 
q  =  0,  and  Eq.  (20).  There  is  excellent  agreement  between  the  analytical  solutions  of  all 
three  curves  and  the  LCS  which  are  found  through  numerical  computation  of  the  FTLE 
flow  fields. 

4.2.  Example  2 — SIS  epidemic  model — External  fluctuations 

As  another  general  application  of  the  extinction  theory  for  finite  populations,  we  con¬ 
sider  the  well-known  problem  of  extinction  in  a  Susceptible-Infectious-Susceptible  (SIS) 
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epidemiological  model.  The  SIS  model  is  given  by  the  following  system  of  equations: 

S  =  /i  —  fiS  +  yl  —  /3IS,  (21a) 

i  =  -(li  +  y)I  +  piS,  (21b) 

where  fi  represents  a  constant  birth  and  death  rate,  /3  represents  the  contact  rate,  and  y 
denotes  the  rate  of  recovery.  If  we  assume  that  the  total  population  size  is  constant  and 
can  be  normalized  to  S  +  I  =  1,  then  Eqs.  (21a)-(21b)  can  be  rewritten  as  the  following 
one-dimensional  (ID)  equation: 

/  = -Ox +  y)/ +  £/(!-/).  (22) 


The  stochastic  version  of  Eq.  (22)  is  given  as 

/  =  -(/*  +  Y)I  +  PI  (I  - 1)  +  ri(t)  =  f  (/)  +  n(t), 


(23) 


where  r\{t)  is  uncorrelated  Gaussian  noise  with  zero  mean  and  models  random  migration 
to  and  from  another  population  (Alonso  et  al.,  2006;  Doering  et  al.,  2005).  Equation  (23) 
has  two  equilibrium  points  given  by  I  =  0  (corresponding  to  the  disease-free  state)  and 
I  =  1  —  (fi  +  y)//3  (corresponding  to  the  endemic  state).  One  can  use  the  Euler-Lagrange 
equation  of  motion  to  find  the  optimal  path  of  extinction  from  the  endemic  state  to  the 
disease-free  state,  where  the  Lagrangian  is  determined  by  Eq.  (23)  and  is  given  as  follows: 

L(I,  /)  =  =  [/  -  (24) 

Computation  of  the  Euler-Lagrange  equation  gives  the  following: 

—  1  =  0.  (25) 


If  one  multiplies  Eq.  (25)  by  I  followed  by  an  integration  with  respect  to  t,  then  one 
obtains 


F(I)2 

2 


(26) 


where  E  is  an  arbitrary  constant.  Using  the  fact  that  the  optimal  path  passes  through  the 
two  equilibrium  points  stated  above,  then  one  finds  that  the  optimal  path  to  extinction  (as 
well  as  its  counterpart  path  from  the  disease-free  state  to  the  endemic  state)  is  given  by 
the  following  equation: 

I  =  ±F(I).  (27) 


As  in  the  first  example,  one  can  numerically  compute  the  optimal  path  to  extinction 
using  the  FTLE.  In  this  example,  we  calculate  the  FTLE  field  using  the  following  2D 
system  which  is  equivalent  to  Eq.  (25): 

i  =  p,  (28a) 

P  =  F(I)F'(I)  =  (pi(  1  -I)-  kI)(P(  1  -  27)  -  k),  (28b) 

where  k  =  i±  +  y .  Figure  3(a)  shows  the  forward  FTLE  plot  computed  using  Eqs.  (28a)- 
(28b)  for  T  =  5,  with  =  5.0  and  k  =  1.0.  In  Fig.  3(a),  one  can  see  that  the  optimal  path 
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Fig.  3  FTLE  flow  fields  computed  using  Eqs.  (28a)-(28b)  with  =  5.0  and  k  =  1.0.  The  integration  time 
is  T  =  5  with  an  integration  step  size  of  t  =  0.1  and  a  grid  resolution  of  0.01  in  both  I  and  p.  (a)  Forward 
FTLE  field  with  the  optimal  path  to  extinction  given  by  the  LCS.  (b)  Average  of  the  forward  and  backward 
FTLE  fields  with  the  optimal  path  to  extinction  and  its  counterpart  optimal  path  given  by  the  LCS  and 
overlaid  with  the  analytical  solution  of  the  optimal  paths  given  by  Eq.  (27).  Note  that  the  averaging  affects 
only  the  value  of  the  FTLE  and  not  the  structure  of  the  FTLE  field.  (Color  figure  online.) 


from  the  endemic  state  to  the  disease-free  state  is  given  by  the  ridge  associated  with  the 
locally  maximal  FTLE. 

One  also  can  find  the  optimal  path  from  the  disease-free  state  to  the  endemic  state  by 
computing  the  backward  FTLE.  By  overlaying  the  forward  and  backward  FTLE  fields, 
one  can  see  the  optimal  path  to  extinction  along  with  its  counterpart  optimal  path  in 
Fig.  3(b).  Also  shown  in  Fig.  3(b)  are  the  two  analytical  solutions  to  the  optimal  path 
to  extinction  and  its  counterpart  optimal  path  which  are  given  by  Eq.  (27).  There  is  excel¬ 
lent  agreement  between  the  analytical  solutions  to  the  two  optimal  paths  and  the  ridges 
which  are  found  through  numerical  computation  of  the  FTLE  flow  fields. 

4.3.  Example  3 — SIS  epidemic  model — Internal  fluctuations 

We  now  consider  the  ID  stochastic  (internal)  version  of  the  SIS  epidemic  model  given  by 
Eq.  (22).  The  probability  Pn(t)  to  observe,  at  time  t,  n  infectious  individuals  is  governed 
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Fig.  4  FTLE  flow  field  computed  using  Eqs.  (31a)-(31b)  with  ft  =  2.0  and  k  =  1.0.  The  integration  time 
is  T  =  10  with  an  integration  step  size  of  t  =  0.1  and  a  grid  resolution  of  0.005  in  both  7  and  p.  The 
optimal  path  to  extinction  is  given  by  the  LCS.  (Color  figure  online.) 


by  the  following  master  equation: 

P,i  =  (M  +  y)[(n  +  l)^n+i  -  nP„]  +  p[(n  -  l)(l  -  (n  -  l))P„_i  -  n(n  -  1  )P„], 

(29) 

Using  the  formalism  of  Gang  (1987),  one  then  has  the  following  Hamiltonian  associated 
with  Eq.  (29): 

H(I,  p)  =  (ji  +  y)l(e~p  -  1)  +  pl{\  -  I)(ep  -  l),  (30) 

and  Hamilton’s  equations  are  therefore  given  as 

/  =  ^-  =  -(ix  +  y)Ie-p  +  0/(1-  I)ep,  (31a) 

dp 

P  =  ~  =  -Qi  +  y)(e~p  -l)+P{ep-  1)(27  -  1).  (31b) 

Although  there  is  no  analytical  solution  for  the  optimal  path  to  extinction  for 
Eqs.  (31a)-(31b),  we  can  once  again  determine  the  optimal  path  by  computing  the  FTLE 
flow  field  associated  with  this  system.  Figure  4  shows  the  forward  FTLE  plot  computed 
using  Eqs.  (31a)-(31b)  for  T  =  10,  with  /3  =  2.0  and  k  =  1.0.  As  we  have  seen  previ¬ 
ously,  the  optimal  path  to  extinction  from  the  endemic  state  to  the  disease-free  state  is 
given  by  the  ridge  associated  with  the  locally  maximal  FTLE. 


5.  Maximal  sensitive  dependence  to  initial  data  near  the  optimal  path 

The  main  heuristic  argument  of  this  section  is  to  show  that  the  path  which  maximizes  the 
probability  of  extinction  also  has  a  finite-time  Lyapunov  exponent  (FTLE)  that  attains  its 
local  maximum  on  the  path.  We  consider  a  general  system  of  equations  and  show  how 
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the  maximum  unstable  direction  along  the  optimal  path  to  extinction  governs  the  local 
hyperbolic  dynamics. 

From  the  Hamiltonian  or  Lagrangian  equations  of  motion,  the  process  which  leads  to 
extinction  consists  of  a  trajectory  which  emanates  from  a  stable  steady  state  xa  g  W1  and 
approaches  the  extinct  state  xs  eR".  Since  the  stable  and  extinct  states  are  both  regular 
saddles  (or  unstable  foci)  in  the  variational  formulations,  they  both  have  hyperbolic  struc¬ 
ture.  Moreover,  every  point  along  the  trajectory  connecting  the  two  states  as  t  =boo  is 
assumed  to  possess  a  local  hyperbolic  structure.  As  an  example,  consider  the  Langevin 
problem  having  a  vector  field  of  position  V  :  W1  —>  Mn,  which  has  the  associated  La¬ 
grangian  L(x ,  x)  =  || jc  —  V(x)||2/2  to  describe  the  action.  Converting  to  a  Hamiltonian 
formulation  leads  one  to  the  following  2n -dimensional  equations  of  motion  and  Hamil¬ 
tonian: 


X  =  p  +  V(x), 

(32a) 

p  =  -V'(x)p, 

(32b) 

\\p\\2 

H(x,p)  =  ^  +  p-V(x), 

(32c) 

where  V'(x)  =  is  the  Jacobian  matrix  evaluated  at  x. 

It  is  immediate  from  Eqs.  (32a)-(32c)  that  {(x,  p)  \  p  =  0}  is  an  invariant  manifold.  In 
addition,  the  optimal  path  must  lie  along  the  H(x,  p)  =  0  surface,  which  means  that  in 
addition  to  the  p  =  0  manifold,  the  zero  surface  includes  {(x,  p)  \  p  =  —2V (x)}. 

To  clarify  the  direction  along  the  optimal  path  as  well  as  the  local  geometry,  we  make 
the  following  assumptions  regarding  V  (x): 

1.  V  (x)  is  smooth, 

2.  V (xa)  =  V (x5)  =  0, 

3.  V'(xfl)  has  eigenvalues  with  negative  real  parts,  and  V'(xs)  has  at  least  one  eigenvalue 
with  positive  real  part. 

Items  2  and  3  imply  that  xa  is  an  attracting  steady  state  and  x>s.  is  an  unstable  steady  state 
in  the  deterministic  dynamical  system.  We  now  assume  that  the  optimal  path  must  satisfy 

lim  (x(t),  p(0)  —  (x„0),  while  lim  (x(t),  p(0)  =  (xfl,  0). 

t — >+oo '  7  t^- oov  7 


Since  H(x(t),  p(t))  =  0  along  the  path,  the  limits  provide  direction  along  the  optimal 
path. 

The  optimal  path  lies  on  the  curve 

C(x,p)  =  |'6  (-00,  oo)  |  pit)  =  -2V(x(0)}, 


and  p  =  0  corresponds  to  the  zero  fluctuation  case.  We  shift  the  optimal  path  to  the  origin 
by  using  the  following  2n -dimensional  transformation: 


u  =  x, 

w  =  p  +  2V(x), 


H(u ,  w)  = 


w 


w  •  V(u). 


(33a) 

(33b) 


2 


(33c) 
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The  new  equations  of  motion  are  now: 


ii  =  dH/dw  =  w  —  V(u), 
w  =  -3H/du  =  V'(u)w. 


(34a) 

(34b) 


The  optimal  path  now  is  described  by  the  curve 

C(U,w)  =  {t  €  (-00,  oo)  |  w(t)  =  0,  u(t)  =  -V(u(t))}, 

while  the  zero  fluctuation  case  given  by  p  =  0  now  corresponds  to  w  =  2V{u). 

The  linearized  variation  along  the  optimal  path  C(U,W)  is  given  by  the  following  matrix 
initial  value  problem  from  Eqs.  (34a)-(34b): 


x=  V'f(r))  X  =  j(u(t),0)X,  X(0  )  =  /. 


(35) 


For  a  fixed  time  to  such  that  u(to)  =  m0,  the  local  eigenvalues  of  Eq.  (35)  are  given  by 
the  eigenvalues  of  ±V'(uo)  and  are  assumed  to  have  nonzero  real  part  for  any  (ii0,  0)  g 
C(u,W)-  Thus,  the  optimal  path  is  hyperbolic  at  every  point.  We  also  suppose  the  existence 
of  a  local  coordinate  system  on  the  path  so  that  there  exists  a  set  of  linearly  independent 
directions  pointwise. 

The  solution  to  the  linear  variational  equation  about  ( u,w)  =  (uq,  0)  for  0  <  t  1  is 
given  by  X(t)  ~  exp  (tJ(uo,  0)).  We  assume  for  simplicity  that  the  eigenvalues  of  V'(mo) 
have  algebraic  multiplicity  of  one.  The  eigenvalues  of  /(i#o,0)  are  given  by  {=LA .,-}”=1, 
where  are  eigenvalues  of  V'(uq).  To  examine  the  dynamic  instability  that  dominates 
locally,  let  A.max  denote  the  eigenvalue  with  largest  real  part  and  be  such  that  Re(\maiX)  >  0. 
Since  —  Amax  has  the  most  negative  real  part,  its  eigendirection  denotes  the  strongest  con¬ 
tracting  direction. 

For  any  (n0,  0)  on  the  path,  the  existence  of  a  set  of  linearly  independent  eigensolutions 
of  J(uo,  0)  implies  there  is  a  transformation  X  =  PY,  with  P  g  L2nx2n  that  diagonalizes 
the  linear  variational  equation  given  by  Eq.  (35).  Without  loss  of  generality,  we  can  also 
assume  that  the  diagonal  consists  of  ordered  descending  eigenvalues  based  on  the  real 
part.  Therefore,  the  linear  variational  system  has  the  form 


■max 


Y  = 


Y. 


(36) 


-max  _ 


For  any  initial  value,  the  solution  to  Eq.  (36)  is 
xp(t;x0)  =  pi  (t),  x2(t),...,  x2„  (0) 


(37) 


=  (e  A10,  e  - A20,  •  •  •  ,  e  XnQ, 

e~X"‘ X(B+1)0,  .  .  •  ,  e~X2t X(2n-l)0,  e~XmmtX2no). 


(38) 
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To  show  that  the  FTLE  takes  it  maximum  along  the  path,  we  notice  that  any  point 
along  the  path  is  hyperbolic  with  a  saddle  structure.  Therefore,  we  consider  an  arbitrary 
initial  condition  lying  within  a  small  domain  containing  the  origin.  Since  almost  any  initial 
condition  hits  the  boundary  of  the  domain  in  finite  time  due  to  the  saddle  structure  of  the 
origin,  we  use  the  escape  time  as  the  final  time  for  the  FTLE.  The  definition  we  use  of  the 
FTLE  is  the  direct  comparison  of  the  distance  between  two  close  trajectories  as  follows: 


cr(t;x  o) 


1 

t 


In  (|  xp(t;x0  +  €) 


xp(t;x  o)  I), 


(39) 


where  e  g  M2n. 

Defining  the  domain  to  be  the  2n -dimensional  hypercube  D  =  [— 1,  l]2n,  then  clearly 
any  point  not  on  the  unstable  manifold  will  escape  in  the  x\  direction  corresponding  to  the 
eigenvalue  with  maximal  real  part.  We  exploit  the  fact  that  the  dynamics  is  governed  by 
the  most  unstable  direction  by  assuming  |A.max|  ^  |,  i  =  2, . . . ,  n.  If  the  initial  condition 

lies  within  a  distance  8  of  the  unstable  manifold  with  0  <  8  1 ,  then  the  time  to  escape 

from  the  domain  for  an  arbitrary  nonzero  initial  condition  is  given  by 


■^-max 


(40) 


Using  the  definition  of  the  exponent  given  by  Eq.  (39),  we  have  that 


O  {t f ,  Xq)  —  ^max  hi 


2  1  k2  1 2 

1  I2 

I  Aa 

+  8  Amax  £2  +  • ' 

■  •  +  5  Amax  €n 

+  5Xmax  €n+i 

I  kl 

H - b  hUmax  £2n_l 


+  \8(:2n  I 


2  In  (8)). 


(41) 


Since  |A.max|  ^  |A;|,  and  since  =bA.max  dominates  the  expanding  and  contracting  direc¬ 
tions,  then  for  8  small,  we  may  just  consider 


v(tf,x  0) 


-Amaxln  (e?/32  +  e|„32) 

2  In  5 


(42) 


Furthermore,  we  find  that 

dajtj^xoi^y)  _  Amaxln(6?)  /  \  /  33  \ 

dS  ~  2S(lnS)2  V  e?  J+  \lnSj' 


(43) 


which  can  be  shown  to  be  negative  assuming  €\  1.  Therefore,  the  FTLE  as  a  function 

of  distance  to  the  stable  invariant  manifold  is  a  decreasing  function,  and  thus  takes  it 
maximum  values  on  the  manifold. 


6.  Conclusions 

In  this  article,  we  have  considered  the  dynamics  of  general  stochastic  epidemic  models 
and  their  extinction  properties  in  finite  populations.  The  random  fluctuations  considered 
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were  from  both  internal  fluctuations,  which  arise  from  mass  action  kinetics,  as  well  as  ex¬ 
ternal  random  forces,  which  may  be  due  to  random  population  migrations.  By  examining 
the  extinction  processes  from  the  master  equation  perspective,  eikonal  approximations  in 
the  large  population  limit  give  a  way  to  solve  for  the  probability  distribution  as  a  func¬ 
tion  of  time.  A  variational  principle  applied  to  the  exponent  of  the  probability  distribution 
near  the  steady  state  was  used  to  maximize  the  probability  to  extinction  from  a  disease 
endemic  state. 

Maximizing  the  probability  to  extinction  means  minimizing  the  action,  which  in  turn 
generates  a  Hamiltonian  (or  Lagrangian)  formulation  that  determines  the  flow  from  en¬ 
demic  to  extinct  states.  The  formulation  describes  the  random  fluctuations  as  a  determin¬ 
istic  effective  force  that  overcomes  the  instability  of  the  extinct  state.  Such  a  deterministic 
flow  describes  the  optimal  path  to  extinction  from  an  endemic  state.  Using  the  above  vari¬ 
ational  formulation,  we  explicitly  derived  the  equations  of  motion  describing  the  optimal 
path  to  extinction  in  three  different  models  from  epidemiology  as  a  two  point  boundary 
value  problem. 

The  main  result  of  our  paper  is  that  the  optimal  path  is  intimately  related  to  the  maxi¬ 
mal  sensitivity  of  the  dynamics  of  unstable  Hamiltonian  flows.  Specifically,  we  introduced 
a  novel  method  to  find  the  optimal  path  to  extinction  that  relies  on  the  computation  of  the 
finite-time  Lyapunov  exponent  field  for  the  dynamical  flow  under  consideration.  The  ex¬ 
ponents  provide  a  measure  of  sensitivity  to  initial  conditions  in  finite  time.  Moreover,  we 
have  shown  that  the  system  possesses  maximal  sensitivity  near  the  optimal  path  to  ex¬ 
tinction.  Therefore,  we  are  able  to  use  the  finite-time  Lyapunov  exponents  to  dynamically 
evolve  toward  the  optimal  path  trajectory. 

To  demonstrate  the  equivalence  of  the  maximal  sensitivity  and  the  optimal  path  max¬ 
imizing  the  probability  of  extinction,  we  have  considered  three  prototypical  examples 
from  mathematical  epidemiology.  In  the  examples,  we  have  considered  both  internal  and 
external  noise,  and  we  have  considered  both  a  Hamiltonian  and  Lagrangian  formulation. 
Furthermore,  in  each  of  the  three  examples,  we  have  shown  that  the  optimal  path  to  ex¬ 
tinction  is  equated  with  having  a  (locally)  maximal  sensitivity  to  initial  condition,  which 
implies  a  relation  at  a  fundamental  level  between  the  optimal  path  and  the  FTLE.  Even 
though  there  exist  many  possible  paths  to  extinction,  the  dynamical  systems  approach 
converges  to  the  path  that  maximizes  the  probability  to  extinction. 

An  example  of  the  evolution  of  the  FTLE  flow  field  showing  the  convergence  of  the 
locally  maximal  FTLE  ridge  to  the  optimal  path  is  shown  in  Fig.  5.  The  FTLE  flow  field 
is  computed  for  T  =  10  using  Eqs.  (31a)-(31b)  for  the  SIS  epidemic  model  with  internal 
fluctuations  (Section  4.3).  Figure  5  shows  snapshots  of  the  FTLE  field  taken  at  t  =  1, 
t  =  3,  t  =  5,  t  =  1,  and  t  =  9.  The  final  snapshot  at  t  =  10  is  shown  in  Fig.  4.  A  video  that 
shows  the  evolution  of  the  FTLE  flow  field  can  be  found  online  as  ancillary  material. 

The  parameter  values  chosen  for  the  three  examples  are  such  that  the  extinct  and  en¬ 
demic  states  are  far  away  from  one  another,  implying  that  the  system  is  operating  far  from 
any  bifurcation  points.  This  result  is  very  important,  since  in  general,  no  approximate 
analytical  treatment,  like  the  one  performed  in  Dykman  et  al.  (2008),  is  possible  if  the 
system’s  dynamics  is  not  sufficiently  close  to  the  bifurcation  point.  However,  any  scaling 
behavior  of  the  exponent  in  the  probability  of  extinction  may  still  be  computed  along  the 
optimal  path,  which  is  the  important  advance  afforded  by  the  new  procedure  proposed  in 
this  paper. 
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Fig.  5  Snapshots  taken  at  (a)  t  =  1,  (b)  t  =  3,  (c)  t  =  5,  (d)  t  =  7,  and  (e)  t  =  9  showing  the  evolution  of 
the  FTLE  flow  field  computed  using  Eqs.  (31a)-(31b)  with  ft  =  2.0  and  k  =  1.0.  The  integration  time  is 
T  =  10  with  an  integration  step  size  of  t  =  0.1  and  a  grid  resolution  of  0.005  in  both  I  and  p.  The  final 
snapshot  in  the  series  is  taken  at  t  =  10  and  is  shown  in  Fig.  4.  A  video  that  shows  the  evolution  of  the 
FTLE  flow  field  can  be  found  online  as  ancillary  material.  (Color  figure  online.) 


In  the  future,  we  plan  on  considering  more  complicated  systems.  Of  particular  interest 
are  bistable  systems  (e.g.,  adaptive  networks  (Shaw  and  Schwartz,  2008)  and  the  Schlogl 
birth-death  process  (Doering  et  al.,  2007)),  and  higher-dimensional  systems  (e.g.,  multi¬ 
strain  epidemic  models  (Shaw  et  al.,  2007)).  Because  the  method  is  general,  and  unifies 
dynamical  systems  theory  with  the  probability  of  extinction,  we  expect  that  any  system 
found  in  other  fields  can  be  understood  using  this  approach. 
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